﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using AdaptiveDeNet.CSV;
using AdaptiveDeNet.RAND;
using AdaptiveDeNet.Functions;

namespace AdaptiveDeNet.DifferentialEvolution
{
    public partial class Optimizer
    {

        private Random systemRandom = new Random();

        /// <summary>
        /// JADE - Adaptive Differential Evolution with external archive,  tries to find global minimum point of given function
        /// </summary>
        /// <param name="lowerBound">Lower bound vector</param>
        /// <param name="upperBound">Upper bound vector</param>
        /// <param name="maxEvals">maximum evaluations of objective function</param>
        /// <param name="dim">dimension of problem</param>
        /// <param name="N">size of population</param>
        /// <param name="myEpps">small nonnegative value, when (Fmax - Fmin) is smaller then myEpps, the termination criteria are satisfied</param>
        /// <param name="c">controls the rate of parameter adaptation, value between 0 and 1</param>
        /// <param name="p">determines the greediness of the mutation strategy (count of selected bestPmembers), value between 0 and 1<</param>
        /// <returns></returns>
        public Result JADE(double[] lowerBound, double[] upperBound, int maxEvals, int dim, int N, double myEpps, double c, double p) 
        {
            //Initial settings

            double uCr = 0.5;
            double uF = 0.5;
            Member[] A = new Member[N];

            this.nRand = new Rand();

            double Fmin = double.MaxValue;
            double Fmax = double.MinValue;

            //deklarace hlavni populace
            Member[] mainPop = new Member[N];
            Member[] trials = new Member[N];

            Member trial;

            //inicializace
            // x = rand(0,1) * (uB - lB) + lb

            double[] uBlBDiff = new double[dim];

            for (int i = 0; i < dim; i++)
            {
                uBlBDiff[i] = upperBound[i] - lowerBound[i];    
            }

            for (int i = 0; i < N; i++)
            {
                Member newMember = new Member(dim);

                // pro kazdou dimenzi
                for (byte j = 0; j < dim; j++)
                {
                    newMember.Vector[j] = nRand.NextDouble() * uBlBDiff[j] + lowerBound[j];
                }

                //evaluace
                newMember.FEval = fEval.Feval(newMember.Vector);
                // zjisteni min. hodnoty  funkce v nulte populaci
                if (newMember.FEval < Fmin) Fmin = newMember.FEval;

                mainPop[i] = newMember;
            }

            Fmax = Common.Common.FindMax(mainPop);
            int popEvals = N;

            //main loop
            while (popEvals < maxEvals && (Fmax - Fmin) > myEpps)
            {
                //set of all successful crossover propabilities
                double[] Scr = new double[N];
                double[] Sf = new double[N];

                 //vygenerovani  
                for (int i = 0; i < N; i++)
                {
                    //for each individual
                    double Cr = MathNet.Numerics.Distributions.Normal.Sample(systemRandom, uCr, 0.1);
                    double F = MathNet.Numerics.Distributions.Cauchy.Sample(systemRandom, uF, 0.1);

                    if (Cr < 0) Cr = 0.0;
                    if (Cr > 1) Cr = 1.0;

                    while (F <= 0) F = MathNet.Numerics.Distributions.Cauchy.Sample(systemRandom, uF, 0.1);
                    if (F >= 1) F = 1;

                    int bestP = GetBestPMemberIndex(mainPop, p);
                    int r1 = this.nRand.Next(N);

                    int PUnionAlength = N + this.PUnionAlength(A);
                    int r2 = this.nRand.Next(PUnionAlength);

                    while (i == r1) r1 = this.nRand.Next(N);
                    while (i == r1 || r1 == r2) r2 = this.nRand.Next(PUnionAlength);

                    int pom = r2;

                    bool liesInAset = false;

                    if (r2 > (N - 1) )
	                {
                        liesInAset = true;
                        r2 = r2 - (N);
	                }

                    trial = new Member(dim);

                    int jrand = this.nRand.Next(dim);

                    //pro kazdou dimenzi
                    for (int j = 0; j < dim; j++)
                    {
                        //krizeni
                        if (this.nRand.NextDouble() <= Cr || j == jrand)
                        {
                            double r2val;

                            if (liesInAset)
                            {
                                r2val = A[r2].Vector[j];
                            }
                            else
                            {
                                r2val = mainPop[r2].Vector[j];
                            }

                            //  
                            // x = Rcurrent + F * (Rpbest - Rcurrent) + F * (r1 - r2)
                            double newVal = mainPop[i].Vector[j] + F * (mainPop[bestP].Vector[j] - mainPop[i].Vector[j])
                                + F * (mainPop[r1].Vector[j] - r2val);


                            if (newVal < lowerBound[j] || newVal > upperBound[j])
                            {
                                newVal = Common.Common.Zrdcad(newVal, lowerBound[j], upperBound[j]);
                                //newVal = nRand.NextDouble() * uBlBDiff[j] + lowerBound[j];
                            }

                            trial.Vector[j] = newVal;
                        }
                        else
                        {
                            trial.Vector[j] = mainPop[i].Vector[j];
                        }
                    }



                    //evaluace  fce
                    trial.FEval = this.fEval.Feval(trial.Vector);
                    trials[i] = trial;
                    popEvals++;

                    if (trials[i].FEval <= mainPop[i].FEval)
                    {
                        this.InsertToA(mainPop[i], A);
                        Scr[i] = Cr;
                        Sf[i] = F;
                    }
                }

                //selekce
                for (int i = 0; i < N; i++)
                {
                    if (trials[i].FEval <= mainPop[i].FEval)
                    {
                        mainPop[i] = trials[i];
                        if (trials[i].FEval < Fmin) Fmin = trials[i].FEval;
                    }
                }

                Fmax = Common.Common.FindMax(mainPop);
                //Fmax = this.FindMax(mainPopulation);
                //Fmin = this.NajdiMin(mainPopulation);

                //uCr - computation
                uCr = (1 - c) * uCr + c * meanA(Scr, popEvals);
                //uF - computation
                uF = (1 - c) * uF + c * meanL(Sf);
            }

            Result result = new Result();

            result.FuncEvals = popEvals;
            result.Fmin = Fmin; 
            result.Point = Common.Common.GetBestMember(mainPop);

            //Console.WriteLine("f : " + Fmin);
            //Console.WriteLine("e : " + popEvals);
            //Console.WriteLine();

            return result;

            


        }

        /// <summary>
        /// JADE - Adaptive Differential Evolution with external archive,  tries to find global minimum point of given function
        /// </summary>
        /// <param name="function">objective function</param>
        /// <param name="lowerBound">Lower bound vector</param>
        /// <param name="upperBound">Upper bound vector</param>
        /// <param name="maxEvals">maximum evaluations of objective function</param>
        /// <param name="dim">dimension of problem</param>
        /// <param name="N">size of population</param>
        /// <param name="myEpps">small nonnegative value, when (Fmax - Fmin) is smaller then myEpps, the termination criteria are satisfied</param>
        /// <param name="c">controls the rate of parameter adaptation, value between 0 and 1</param>
        /// <param name="p">determines the greediness of the mutation strategy (count of selected bestPmembers), value between 0 and 1<</param>
        /// <returns></returns>
        public Result JADE(Func<double[], double> function, double[] lowerBound, double[] upperBound, int maxEvals, int dim, int N, double myEpps, double c, double p)
        {
            //Initial settings

            double uCr = 0.5;
            double uF = 0.5;
            Member[] A = new Member[N];

            this.nRand = new Rand();

            double Fmin = double.MaxValue;
            double Fmax = double.MinValue;

            //deklarace hlavni populace
            Member[] mainPop = new Member[N];
            Member[] trials = new Member[N];

            Member trial;

            //inicializace
            // x = rand(0,1) * (uB - lB) + lb

            double[] uBlBDiff = new double[dim];

            for (int i = 0; i < dim; i++)
            {
                uBlBDiff[i] = upperBound[i] - lowerBound[i];
            }

            for (int i = 0; i < N; i++)
            {
                Member newMember = new Member(dim);

                // pro kazdou dimenzi
                for (byte j = 0; j < dim; j++)
                {
                    newMember.Vector[j] = nRand.NextDouble() * uBlBDiff[j] + lowerBound[j];
                }

                //evaluace
                newMember.FEval = function.Invoke(newMember.Vector);
                // zjisteni min. hodnoty  funkce v nulte populaci
                if (newMember.FEval < Fmin) Fmin = newMember.FEval;

                mainPop[i] = newMember;
            }

            Fmax = Common.Common.FindMax(mainPop);
            int popEvals = N;

            //main loop
            while (popEvals < maxEvals && (Fmax - Fmin) > myEpps)
            {
                //set of all successful crossover propabilities
                double[] Scr = new double[N];
                double[] Sf = new double[N];

                //vygenerovani  
                for (int i = 0; i < N; i++)
                {
                    //for each individual
                    double Cr = MathNet.Numerics.Distributions.Normal.Sample(systemRandom, uCr, 0.1);
                    double F = MathNet.Numerics.Distributions.Cauchy.Sample(systemRandom, uF, 0.1);

                    if (Cr < 0) Cr = 0.0;
                    if (Cr > 1) Cr = 1.0;

                    while (F <= 0) F = MathNet.Numerics.Distributions.Cauchy.Sample(systemRandom, uF, 0.1);
                    if (F >= 1) F = 1;

                    int bestP = GetBestPMemberIndex(mainPop, p);
                    int r1 = this.nRand.Next(N);

                    int PUnionAlength = N + this.PUnionAlength(A);
                    int r2 = this.nRand.Next(PUnionAlength);

                    while (i == r1) r1 = this.nRand.Next(N);
                    while (i == r1 || r1 == r2) r2 = this.nRand.Next(PUnionAlength);

                    int pom = r2;

                    bool liesInAset = false;

                    if (r2 > (N - 1))
                    {
                        liesInAset = true;
                        r2 = r2 - (N);
                    }

                    trial = new Member(dim);

                    int jrand = this.nRand.Next(dim);

                    //pro kazdou dimenzi
                    for (int j = 0; j < dim; j++)
                    {
                        //krizeni
                        if (this.nRand.NextDouble() <= Cr || j == jrand)
                        {
                            double r2val;

                            if (liesInAset)
                            {
                                r2val = A[r2].Vector[j];
                            }
                            else
                            {
                                r2val = mainPop[r2].Vector[j];
                            }

                            //  
                            // x = Rcurrent + F * (Rpbest - Rcurrent) + F * (r1 - r2)
                            double newVal = mainPop[i].Vector[j] + F * (mainPop[bestP].Vector[j] - mainPop[i].Vector[j])
                                + F * (mainPop[r1].Vector[j] - r2val);


                            if (newVal < lowerBound[j] || newVal > upperBound[j])
                            {
                                newVal = Common.Common.Zrdcad(newVal, lowerBound[j], upperBound[j]);
                                //newVal = nRand.NextDouble() * uBlBDiff[j] + lowerBound[j];
                            }

                            trial.Vector[j] = newVal;
                        }
                        else
                        {
                            trial.Vector[j] = mainPop[i].Vector[j];
                        }
                    }



                    //evaluace  fce
                    trial.FEval = function.Invoke(trial.Vector);
                    trials[i] = trial;
                    popEvals++;

                    if (trials[i].FEval <= mainPop[i].FEval)
                    {
                        this.InsertToA(mainPop[i], A);
                        Scr[i] = Cr;
                        Sf[i] = F;
                    }
                }

                //selekce
                for (int i = 0; i < N; i++)
                {
                    if (trials[i].FEval <= mainPop[i].FEval)
                    {
                        mainPop[i] = trials[i];
                        if (trials[i].FEval < Fmin) Fmin = trials[i].FEval;
                    }
                }

                Fmax = Common.Common.FindMax(mainPop);
                //Fmax = this.FindMax(mainPopulation);
                //Fmin = this.NajdiMin(mainPopulation);

                //uCr - computation
                uCr = (1 - c) * uCr + c * meanA(Scr, popEvals);
                //uF - computation
                uF = (1 - c) * uF + c * meanL(Sf);
            }

            Result result = new Result();

            result.FuncEvals = popEvals;
            result.Fmin = Fmin;
            result.Point = Common.Common.GetBestMember(mainPop);

            //Console.WriteLine("f : " + Fmin);
            //Console.WriteLine("e : " + popEvals);
            //Console.WriteLine();

            return result;




        }

        public CEC14Result JADE(double[] lowerBound, double[] upperBound, int maxEvals, int dim, int N, double myEpps, double c, double p, Queue<double> FESMultipliers, double fX, int funcNum)
        {
            CEC14Result result = new CEC14Result();

            //Initial settings

            double uCr = 0.5;
            double uF = 0.5;
            Member[] A = new Member[N];

            this.nRand = new Rand();

            double Fmin = double.MaxValue;
            double Fmax = double.MinValue;

            //deklarace hlavni populace
            Member[] mainPop = new Member[N];
            Member[] trials = new Member[N];

            Member trial;

            //inicializace
            // x = rand(0,1) * (uB - lB) + lb

            double[] uBlBDiff = new double[dim];

            for (int i = 0; i < dim; i++)
            {
                uBlBDiff[i] = upperBound[i] - lowerBound[i];
            }

            for (int i = 0; i < N; i++)
            {
                Member newMember = new Member(dim);

                // pro kazdou dimenzi
                for (byte j = 0; j < dim; j++)
                {
                    newMember.Vector[j] = nRand.NextDouble() * uBlBDiff[j] + lowerBound[j];
                }

                //evaluace
                newMember.FEval = fEval.Feval(newMember.Vector);
                // zjisteni min. hodnoty  funkce v nulte populaci
                if (newMember.FEval < Fmin) Fmin = newMember.FEval;

                mainPop[i] = newMember;
            }

            Fmax = Common.Common.FindMax(mainPop);
            int popEvals = N;

            //main loop
            while (popEvals < maxEvals)
            {
                //set of all successful crossover propabilities
                double[] Scr = new double[N];
                double[] Sf = new double[N];

                //vygenerovani  
                for (int i = 0; i < N; i++)
                {
                    //for each individual
                    double Cr = MathNet.Numerics.Distributions.Normal.Sample(systemRandom, uCr, 0.1);
                    double F = MathNet.Numerics.Distributions.Cauchy.Sample(systemRandom, uF, 0.1);

                    if (Cr < 0) Cr = 0.0;
                    if (Cr > 1) Cr = 1.0;

                    while (F <= 0) F = MathNet.Numerics.Distributions.Cauchy.Sample(systemRandom, uF, 0.1);
                    if (F >= 1) F = 1;

                    int bestP = GetBestPMemberIndex(mainPop, p);
                    int r1 = this.nRand.Next(N);

                    int PUnionAlength = N + this.PUnionAlength(A);
                    int r2 = this.nRand.Next(PUnionAlength);

                    while (i == r1) r1 = this.nRand.Next(N);
                    while (i == r1 || r1 == r2) r2 = this.nRand.Next(PUnionAlength);

                    int pom = r2;

                    bool liesInAset = false;

                    if (r2 > (N - 1))
                    {
                        liesInAset = true;
                        r2 = r2 - (N);
                    }

                    trial = new Member(dim);

                    int jrand = this.nRand.Next(dim);

                    //pro kazdou dimenzi
                    for (int j = 0; j < dim; j++)
                    {
                        //krizeni
                        if (this.nRand.NextDouble() <= Cr || j == jrand)
                        {
                            double r2val;

                            if (liesInAset)
                            {
                                r2val = A[r2].Vector[j];
                            }
                            else
                            {
                                r2val = mainPop[r2].Vector[j];
                            }

                            //  
                            // x = Rcurrent + F * (Rpbest - Rcurrent) + F * (r1 - r2)
                            double newVal = mainPop[i].Vector[j] + F * (mainPop[bestP].Vector[j] - mainPop[i].Vector[j])
                                + F * (mainPop[r1].Vector[j] - r2val);


                            if (newVal < lowerBound[j] || newVal > upperBound[j])
                            {
                                newVal = Common.Common.Zrdcad(newVal, lowerBound[j], upperBound[j]);
                                //newVal = nRand.NextDouble() * uBlBDiff[j] + lowerBound[j];
                            }

                            trial.Vector[j] = newVal;
                        }
                        else
                        {
                            trial.Vector[j] = mainPop[i].Vector[j];
                        }
                    }



                    //evaluace  fce
                    trial.FEval = this.fEval.Feval(trial.Vector);
                    trials[i] = trial;
                    popEvals++;

                    if (trials[i].FEval <= mainPop[i].FEval)
                    {
                        this.InsertToA(mainPop[i], A);
                        Scr[i] = Cr;
                        Sf[i] = F;
                    }
                }

                //selekce
                for (int i = 0; i < N; i++)
                {
                    if (trials[i].FEval <= mainPop[i].FEval)
                    {
                        mainPop[i] = trials[i];
                        if (trials[i].FEval < Fmin) Fmin = trials[i].FEval;
                    }
                }

                Fmax = Common.Common.FindMax(mainPop);
                //Fmax = this.FindMax(mainPopulation);
                //Fmin = this.NajdiMin(mainPopulation);

                //uCr - computation
                uCr = (1 - c) * uCr + c * meanA(Scr, popEvals);
                //uF - computation
                uF = (1 - c) * uF + c * meanL(Sf);

                if (popEvals == (FESMultipliers.Peek() * maxEvals))
                {
                    result.Records.Add(Fmin - fX);
                    FESMultipliers.Dequeue();
                }
            }


            result.FuncEvals = popEvals;

            //Console.WriteLine("f : " + Fmin);
            //Console.WriteLine("e : " + popEvals);
            //Console.WriteLine();

            return result;
        }

        private void InsertToA(Member member, Member[] A)
        {
            bool wasInserted = false;

            //insert on first null occurence
            for (int i = 0; i < A.Length; i++)
            {
                if (A[i] == null)
                {
                    A[i] = member;
                    wasInserted = true;
                }    
            }

            if (!wasInserted)
            {
                A[this.nRand.Next(A.Length)] = member;    
            }
        }


        private int PUnionAlength(Member[] A)
        {
            int len = 0;

            for (int i = 0; i < A.Length; i++)
            {
                if (A[i] != null)
                {
                    len++;
                }
            }

            return len;
        }

        private int GetBestPMemberIndex(Member[] mainPop, double p)
        {
            int count = (int)Math.Round(mainPop.Length * p);

            int[] bestPMembersIndicies = new int[count];

            double latestMin = double.MinValue;

            for (int i = 0; i < count; i++)
            {
                double minVal = double.MaxValue;

                for (int j = 0; j < mainPop.Length; j++)
                {
                    if (minVal > mainPop[j].FEval && mainPop[j].FEval > latestMin)
                    {
                        bestPMembersIndicies[i] = j;
                        minVal = mainPop[j].FEval;
                    }
                }

                latestMin = minVal;
            }

            return bestPMembersIndicies[this.nRand.Next(count)];
        }


        private double meanL(double[] Sf)
        {
            double f2 = 0;

            for (int i = 0; i < Sf.Length; i++)
            {
                f2 += (Sf[i] * Sf[i]);    
            }

            double f = 0;

            for (int i = 0; i < Sf.Length; i++)
            {
                f += Sf[i];
            }

            double ret = f2 / f;

            if (double.IsNaN(ret) || double.IsInfinity(ret))
            {
                return 0;
            }
            else
            {
                return ret;
            }
         
        }

        private double meanA(double[] Scr, int popEvals)
        {
            double mean = 0;
            double len = 0;

            for (int i = 0; i < Scr.Length; i++)
            {
                mean += Scr[i];

                if (Scr[i] != 0)
                {
                    len++;    
                }
            }

            if (len == 0 || mean == 0)
            {
                return 0;   
            }

            return mean / len;
        }

    }
}
